#! /usr/bin/python3
# -*- coding: UTF-8 -*-
  
'''※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※
File Name: arithmetic_analysis.py
Author: GID5564
Description: 方程求解，算术分析
Version: 1.0
Created Time: 26/11/24-13:49:15
※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※※'''

import math
import numpy as np
from numpy import array, cos, sin, radians, cross  # type: ignore
from typing import List

from sympy import diff
from decimal import Decimal, getcontext

def ucal(u, p):
    """
    计算给定参数的特定数学表达式。
    
    参数:
    u (float): 基数，用于计算的基础值。
    p (int): 幂次，表示乘法操作的次数。
    
    返回:
    float: 计算结果，根据公式 u * (u-1) * (u-2) * ... * (u-(p-1)) 得到的结果。
    
    示例:
    >>> ucal(1, 2)
    0
    >>> ucal(1.1, 2)
    0.11000000000000011
    >>> ucal(1.2, 2)
    0.23999999999999994
    """
    temp = u  # 初始化临时变量 temp 为基数 u
    for i in range(1, p):  # 从 1 到 p-1 进行循环
        temp = temp * (u - i)  # 每次循环将 temp 乘以 (u - i)
    return temp  # 返回最终计算结果

def newton(function, function1, startingInt):
    """
    使用牛顿法求解方程的根。

    参数:
    function (callable): 目标函数，即需要求其根的函数。
    function1 (callable): 目标函数的导数。
    startingInt (float): 迭代的初始值。

    返回:
    float: 方程的近似根。
    """
    x_n = startingInt  # 初始化迭代值
    while True:
        # 计算下一个迭代值
        x_n1 = x_n - function(x_n) / function1(x_n)
        # 如果新旧迭代值之差小于设定的阈值，则认为已收敛，返回结果
        if abs(x_n - x_n1) < 10 ** -5:
            return x_n1
        # 更新迭代值
        x_n = x_n1
        
"""
#例:
def f(x):
    return (x ** 3) - (2 * x) - 5

def f1(x):
    return 3 * (x ** 2) - 2

if __name__ == "__main__":
    print(newton(f, f1, 3))
"""


def intersection(function, x0, x1):
    """
    使用割线法（Secant Method）求解函数的根。

    参数:
    function (callable): 需要找到其根的函数。
    x0 (float): 第一个初始猜测值。
    x1 (float): 第二个初始猜测值。

    返回:
    float: 函数的近似根。
    """
    # 初始化变量，x_n 和 x_n1 分别为两个初始猜测值
    x_n = x0
    x_n1 = x1
    
    while True:
        # 计算新的近似根 x_n2
        x_n2 = x_n1 - (function(x_n1) / ((function(x_n1) - function(x_n)) / (x_n1 - x_n)))
        
        # 如果新旧近似根之间的差值小于设定的阈值，则认为已经收敛，返回结果
        if abs(x_n2 - x_n1) < 10 ** -5:
            return x_n2
        
        # 更新 x_n 和 x_n1 为新的近似根，准备下一次迭代
        x_n = x_n1
        x_n1 = x_n2

"""
#例:
def f(x):
    return math.pow(x, 3) - (2 * x) - 5
if __name__ == "__main__":
    print(intersection(f, 3, 3.5))

"""


def bisection(function, a, b):
    """
    使用二分法在区间 [a, b] 内找到函数的根。
    
    参数:
    function (callable): 需要找到根的函数。
    a (float): 区间的左端点。
    b (float): 区间的右端点。
    
    返回值:
    float: 函数在区间 [a, b] 内的根，如果找不到则返回 None。
    """
    
    start = a  # 初始化起始点为 a
    end = b  # 初始化结束点为 b
    
    if function(a) == 0:  # 如果 a 是函数的根
        return a
    elif function(b) == 0:  # 如果 b 是函数的根
        return b
    elif function(a) * function(b) > 0:  # 如果 a 和 b 同号，则无法找到根
        print("couldn't find root in [a,b]")
        return None
    else:
        mid = start + (end - start) / 2.0  # 计算中点
        while abs(start - mid) > 10 ** -7:  # 当区间长度大于 10^-7 时继续迭代
            if function(mid) == 0:  # 如果中点是根
                return mid
            elif function(mid) * function(start) < 0:  # 如果中点和起点异号
                end = mid  # 更新结束点为中点
            else:  # 如果中点和起点同号
                start = mid  # 更新起始点为中点
            mid = start + (end - start) / 2.0  # 重新计算中点
        return mid  # 返回最终的近似根

"""
#例:
def f(x):
    return math.pow(x, 3) - 2 * x - 5

if __name__ == "__main__":
    print(bisection(f, 1, 1000))
"""


def retroactive_resolution(coefficients: np.matrix, vector: np.array) -> np.array:
    """
    该函数执行一个用于三角矩阵的后向线性系统求解。

    示例：
        2x1 + 2x2 - 1x3 = 5         2x1 + 2x2 = -1
        0x1 - 2x2 - 1x3 = -7        0x1 - 2x2 = -1
        0x1 + 0x2 + 5x3 = 15
    >>> retroactive_resolution(np.matrix([[2, 2, -1], [0, -2, -1], [0, 0, 5]]), np.array([5, -7, 15]))
    array([[2.],
           [2.],
           [3.]])
    >>> retroactive_resolution(np.matrix([[2, 2], [0, -2]]), np.array([-1, -1]))
    array([[-1. ],
           [ 0.5]])
    """

    # 获取系数矩阵的行数和列数
    rows, columns = np.shape(coefficients)

    # 初始化解向量，全为0
    x = np.zeros((rows, 1), dtype=float)
    
    # 从最后一行开始向上迭代
    for row in reversed(range(rows)):
        sum = 0
        # 计算当前行的已知变量之和
        for col in range(row + 1, columns):
            sum += coefficients[row, col] * x[col]

        # 更新当前行的解
        x[row, 0] = (vector[row] - sum) / coefficients[row, row]

    return x


def gaussian_elimination(coefficients: np.matrix, vector: np.array) -> np.array:
    """
    使用高斯消元法解线性方程组。

    参数:
    coefficients (np.matrix): 系数矩阵，必须为方阵。
    vector (np.array): 常数项向量。

    返回值:
    np.array: 方程组的解向量。

    示例:
    >>> gaussian_elimination([[1, -4, -2], [5, 2, -2], [1, -1, 0]], [[-2], [-3], [4]])
    array([[ 2.3 ],
           [-1.7 ],
           [ 5.55]])
    >>> gaussian_elimination([[1, 2], [5, 2]], [[5], [5]])
    array([[0. ],
           [2.5]])
    """
    # 检查输入的系数矩阵是否为方阵
    rows, columns = np.shape(coefficients)
    if rows != columns:
        raise ValueError("系数矩阵必须是方阵")

    # 构建增广矩阵
    augmented_mat = np.concatenate((coefficients, vector), axis=1)
    augmented_mat = augmented_mat.astype("float64")

    # 将矩阵转换为上三角矩阵
    for row in range(rows - 1):
        pivot = augmented_mat[row, row]
        if pivot == 0:
            raise ValueError("主元不能为零")
        for col in range(row + 1, columns):
            factor = augmented_mat[col, row] / pivot
            augmented_mat[col, :] -= factor * augmented_mat[row, :]

    # 回代求解
    x = retroactive_resolution(
        augmented_mat[:, 0:columns], augmented_mat[:, columns : columns + 1]
    )

    return x
    
    
    

def polar_force(
    magnitude: float, angle: float, radian_mode: bool = False
) -> List[float]:
    """
    将极坐标下的力分解为直角坐标系下的分量。
    参数:
        magnitude (float): 力的幅值（大小）。
        angle (float): 力的角度，单位可以是度或弧度。
        radian_mode (bool): 如果为True，则角度以弧度为单位；否则以度为单位。
    返回:
        List[float]: 力的x和y分量。
    示例:
        >>> polar_force(10, 45)
        [7.0710678118654755, 7.071067811865475]
        >>> polar_force(10, 3.14, radian_mode=True)
        [-9.999987317275394, 0.01592652916486828]
    """
    if radian_mode:
        # 如果角度是以弧度表示的，直接使用cos和sin函数计算x和y分量
        return [magnitude * cos(angle), magnitude * sin(angle)]
    # 如果角度是以度表示的，先转换为弧度再计算x和y分量
    return [magnitude * cos(radians(angle)), magnitude * sin(radians(angle))]

def in_static_equilibrium(
    forces: np.ndarray, location: np.ndarray, eps: float = 10 ** -1
) -> bool:
    """
    检查一个系统是否处于平衡状态。
    参数:
        forces (np.ndarray): 力的数组，每个元素是一个包含两个元素的列表，分别表示力的x和y分量。
        location (np.ndarray): 位置的数组，每个元素是一个包含两个元素的列表，分别表示位置的x和y坐标。
        eps (float): 判断平衡状态的容差。
    返回:
        bool: 如果系统处于平衡状态，返回True；否则返回False。
    示例:
        >>> force = np.array([[1, 1], [-1, 2]])
        >>> location = np.array([[1, 0], [10, 0]])
        >>> in_static_equilibrium(force, location)
        False
    """
    # 计算力矩的和，力矩是位置向量和力向量的叉积
    moments: np.ndarray = np.cross(location, forces)
    sum_moments: float = np.sum(moments)
    # 如果力矩的和小于给定的容差，则认为系统处于平衡状态
    return abs(sum_moments) < eps

"""
#例:
if __name__ == "__main__":

    # Test to check if it works

    forces = array(
        [polar_force(718.4, 180 - 30), polar_force(879.54, 45), polar_force(100, -90)]
    )

    location = array([[0, 0], [0, 0], [0, 0]])

    assert in_static_equilibrium(forces, location)

    # Problem 1 in image_data/2D_problems.jpg
    forces = array(
        [
            polar_force(30 * 9.81, 15),
            polar_force(215, 180 - 45),
            polar_force(264, 90 - 30),
        ]
    )

    location = array([[0, 0], [0, 0], [0, 0]])

    assert in_static_equilibrium(forces, location)

    # Problem in image_data/2D_problems_1.jpg
    forces = array([[0, -2000], [0, -1200], [0, 15600], [0, -12400]])

    location = array([[0, 0], [6, 0], [10, 0], [12, 0]])

    assert in_static_equilibrium(forces, location)

    import doctest

    doctest.testmod()
"""

    

def LUDecompose(table):
    """
    对给定的方阵进行LU分解。Lower-Upper (LU) Decomposition.
    
    参数:
    table (numpy.ndarray): 需要进行LU分解的方阵。
    
    返回值:
    tuple: 包含两个元素的元组，第一个元素是下三角矩阵L，第二个元素是上三角矩阵U。
           如果输入的矩阵不是方阵，则返回空列表。
    """
    # 获取矩阵的行数和列数
    rows, columns = numpy.shape(table)
    
    # 初始化L和U矩阵为零矩阵
    L = numpy.zeros((rows, columns))
    U = numpy.zeros((rows, columns))
    
    # 检查输入是否为方阵
    if rows != columns:
        return []
    
    # 开始LU分解过程
    for i in range(columns):
        # 计算L矩阵的第i行第j列的元素
        for j in range(i):
            sum = 0
            for k in range(j):
                sum += L[i][k] * U[k][j]
            L[i][j] = (table[i][j] - sum) / U[j][j]
        
        # 设置L矩阵的对角线元素为1
        L[i][i] = 1
        
        # 计算U矩阵的第i行第j列的元素
        for j in range(i, columns):
            sum1 = 0
            for k in range(i):
                sum1 += L[i][k] * U[k][j]
            U[i][j] = table[i][j] - sum1
    
    # 返回L和U矩阵
    return L, U

"""
#例:
if __name__ == "__main__":
    matrix = numpy.array([[2, -2, 1], [0, 1, 2], [5, 3, 1]])
    L, U = LUDecompose(matrix)
    print(L)
    print(U)
"""
    


def newton_raphson(func: str, a: int, precision: int = 10 ** -10) -> float:
    """
    使用牛顿-拉夫森方法从点 'a' 开始寻找函数的根。
    
    参数:
    func (str): 需要求根的函数，以字符串形式表示。
    a (int): 初始猜测值。
    precision (float): 精度要求，默认为 10^-10。
    
    返回:
    float: 找到的根。
    
    示例:
    >>> newton_raphson("sin(x)", 2)
    3.1415926536808043
    >>> newton_raphson("x**2 - 5*x +2", 0.4)
    0.4384471871911695
    >>> newton_raphson("x**2 - 5", 0.1)
    2.23606797749979
    >>> newton_raphson("log(x)- 1", 2)
    2.718281828458938
    """
    # 设置计算精度
    getcontext().prec = 28
    
    # 将初始猜测值转换为Decimal类型
    x = Decimal(a)
    
    while True:
        # 计算当前点的函数值和导数值
        f_val = eval(func.replace('x', f'({x})'))
        f_prime_val = eval(f"{func.replace('x', '')}_prime".replace('x', f'({x})'))
        
        # 更新x的值
        x = x - Decimal(f_val) / Decimal(f_prime_val)
        
        # 如果当前值与上一次迭代值的差小于精度要求，则认为已经找到根
        if abs(Decimal(f_val)) < precision:
            return float(x)

"""
#倒:
if __name__ == "__main__":
    # Find root of trigonometric function
    # Find value of pi
    print(f"The root of sin(x) = 0 is {newton_raphson('sin(x)', 2)}")
    # Find root of polynomial
    print(f"The root of x**2 - 5*x + 2 = 0 is {newton_raphson('x**2 - 5*x + 2', 0.4)}")
    # Find Square Root of 5
    print(f"The root of log(x) - 1 = 0 is {newton_raphson('log(x) - 1', 2)}")
    # Exponential Roots
    print(f"The root of exp(x) - 1 = 0 is {newton_raphson('exp(x) - 1', 0)}")

"""
    
    
def SecantMethod(lower_bound, upper_bound, repeats):
    """
    使用割线法求解方程的根。
    
    参数:
    lower_bound (float): 区间的下界。
    upper_bound (float): 区间的上界。
    repeats (int): 迭代次数。
    
    返回:
    float: 近似的根值。
    
    示例:
    >>> SecantMethod(1, 3, 2)
    0.2139409276214589
    """
    x0 = lower_bound  # 初始化第一个近似值
    x1 = upper_bound  # 初始化第二个近似值
    for i in range(repeats):  # 进行指定次数的迭代
        # 更新近似值，使用割线法公式
        x0, x1 = x1, x1 - (f(x1) * (x1 - x0)) / (f(x1) - f(x0))
    return x1  # 返回最终的近似根值

"""
#例:
def f(x):

    """
    >>> f(5)
    39.98652410600183
    """
    return 8 * x - 2 * exp(-x)
    
print(f"The solution is: {SecantMethod(1, 3, 2)}")

"""
    
    
    
    
    
    
    
    
    
    
    




















